arXiv: 1503.05914v2 [physics.soc-ph] 8 Apr 2016 


Macroscopic description of complex adaptive networks co-evolving with dynamic node 

states 

Marc Wiedermann,^’^’[^ Jonathan F. Donges/’^ Jobst Heitzig,^ Wolfgang Lucht/’^ and Jurgen Kurths^’^’^’^ 

^Potsdam Institute for Climate Impaet Researeh — P.O. Box 60 12 03, 144^^ Potsdam, Germany, EU 
^Department of Physics, Humboldt University — Newtonstr. 15, 12489 Berlin, Germany, EU 
^Stockholm Resilience Centre, Stockholm University — Kraftriket 2B, II 4 19 Stockholm, Sweden, EU 
^Department of Geography, Humboldt University — Rudower Chaussee 16, 12489 Berlin, Germany, EU 
^Institute for Complex Systems and Mathematical Biology, 

University of Aberdeen — Aberdeen AB24 SEX, UK, EU 

^Department of Control Theory, Nizhny Novgorod State University — Gagarin Avenue 23, 606950 Nizhny Novgorod, Russia 

(Dated: April 11, 2016) 

In many real-world complex systems, the time-evolution of the network’s structure and the dy¬ 
namic state of its nodes are closely entangled. Here, we study opinion formation and imitation 
on an adaptive complex network which is dependent on the individual dynamic state of each node 
and vice versa to model the co-evolution of renewable resources with the dynamics of harvesting 
agents on a social network. The adaptive voter model is coupled to a set of identical logistic growth 
models and we mainly find that in such systems, the rate of interactions between nodes as well 
as the adaptive rewiring probability are crucial parameters for controlling the sustainability of the 
system’s equilibrium state. We derive a macroscopic description of the system which provides a 
general framework to model and quantify the influence of single node dynamics on the macroscopic 
state of the network. The thus obtained framework is applicable to many fields of study, such as 
epidemic spreading, opinion formation or socio-ecological modeling. 

PACS numbers: 89.75.Fb, 89.75.Hc, 89.65.-s, 87.23.Ge 


I. INTRODUCTION 

Complex network theory has proven to be a powerful 
tool for studying properties, dynamics and evolution of 
many real-world complex systems m- Of particular in¬ 
terest is to investigate adaptive or temporal networks and 
their respective dynamics laHS- Typical processes stud¬ 
ied in this field are epidemic spreading [MSI or opinion 
formation, e.g., based on the adaptive voter model mm- 
Interactions are modeled by randomly picking a pair of 
linked nodes and, with fixed probabilities, either chang¬ 
ing the state of one of the two nodes or modifying their 
neighborhood structure by adaptive rewiring. However, 
recent results have emphasized that opinion formation 
and imitation processes in fact do not take place with 
fixed probabilities but can depend on the payoff or per¬ 
formance of different opinion-related choices made by the 
agents or nodes involved [IIHI3]. 

In addition to the structure and dynamics of networks 
there has been a variety of studies on the dynamics on 
networks, where nodes in the network represent indi¬ 
vidual dynamical systems and links indicate directed or 
symmetric interactions between them [H US- It has 
been suggested that the interplay between the dynamics 
of and on networks should be much more investigated, 
since the dynamics of each of the coupled subsystems is 
expected to change significantly when compared to their 
autonomous time-evolution [3]. 
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In this work, we propose a model that combines both 
aspects. For this purpose we refine the adaptive voter 
model so that there is no fixed probability for pairs of 
nodes to either imitate each other’s opinion or adaptively 
rewire their acquaintance structure. Instead, each node 
also represents a dynamical system which, for illustra¬ 
tion, is chosen here to be simple and easily understood 
if treated in an isolated fashion. In particular, we choose 
a logistic growth model, which is a paradigm for the dy¬ 
namics of a bounded renewable resource m- Whenever 
interactions between nodes take place, the states of the 
respective dynamical systems are also taken into account. 
As a consequence, imitation processes depend explicitly 
on the nodes’ states as well as on the current network 
structure. At the same time each of the nodes’ opinions 
influences a parameter of the local dynamical system. 

The proposed model serves as a narrative for possi¬ 
bly emerging dynamics in co-evolutionary human-nature 
interactions [MIS- It complements conceptual stud¬ 
ies on the effects of economic growth on the ecospheric 
state [ 20 I Hi] as well as work on resource exploitation 
models that take into account the co-evolution of stylized 
resource dynamics with a similarly paradigmatic popu¬ 
lation growth model [22l|23]. The proposed model, for 
the first time, takes into account individual pairwise in¬ 
teractions of agents on a social network when studying 
the stability and dynamics of such intertwined systems. 

So far, in the context of sustainability science [24] . 
studies on the effect of different exploitation strategies 
on the state of a certain ecospheric component have been 
carried out by, e.g., studying the extraction of water in 
rivers by a network of interconnected harvesters [25H27| . 
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However, no systematic analysis of the underlying net¬ 
work structure and resulting dynamics was performed. 
In addition, no network dynamics, such as adaptation or 
imitation processes, were included in these studies and 
the focus was mainly set on studying the state of the eco- 
sphere for different harvesting strategies that were evolv¬ 
ing deterministically in order to optimize all harvesters’ 
payoffs. 

In contrast, imitation dynamics with high numbers of 
agents or players have been studied in the context of 
evolutionary game theory [121 [ISl UHl [22] • However, in 
no such cases the dynamics of resources or other exter¬ 
nalities has been taken into account and, hence, no co¬ 
evolution of different subsystems has been studied. Here, 
the proposed model serves to illustrate the rich dynam¬ 
ics that may emerge from the coupling of these different 
subsystems, even though the complexity in each of the 
subcomponents remains manageable. 

After the introduction of all key components and pro¬ 
cesses constituting the model in Sec. [n| we perform nu¬ 
merical simulations of the system. In Sec. in we first 
study the case of a static network and no is adaptation 
taking place. We find that the system converges into 
either a state where all logistic growth models, e.g., re¬ 
sources, converge into a state of full depletion or into a 
state of positive stock. The latter is to be interpreted 
as the more sustainable and, hence, desired outcome of 
the model. We uncover that the likelihood to converge 
into either of the two states is mainly determined by the 
frequency of interactions between nodes. 

In Sec. [IV| we then study the effect of network adap¬ 
tation and show that the stability of the system changes 
in dependence on the choice of the adaptation frequency. 
Specifically we deduct that for each interaction frequency 
there exists an appropriate rate of network adaptation, 
such that the system can be guided into a sustainable 
state. 

Finally, we derive a low-dimensional set of rate equa¬ 
tions for variables that approximate the model’s macro¬ 
scopic state in Sec. HIB for the static and in Sec. |IV| for 
the adaptive case. These equations are generally appli¬ 
cable to any study of opinion formation or spreading if 
the probabilities of changes in node states by imitation 
are appropriately chosen. Finally, conclusion are drawn 
in Sec. El 


II. MODEL DESCRIPTION 

Assume a temporal network G{V,L{t)) consisting of a 
fixed set of N nodes V = {vi,V 2 ,..., and an evolving 
set of links L{t). It is represented by the time-dependent 
adjacency matrix A{t). Each node Vi represents a re¬ 
newable resource stock Si{t) that obeys a logistic growth 
model and is harvested with an effort level Ei{t) [16], 

= aiSi{t){l - Si{t)/Ki) - qiSi{t)Ei{t). ( 1 ) 


For this study, we set the growth rates = 1, capacities 
Ki = 1 and catch-coefficients = 1 for alH = 1,..., 
and measure the time and stocks in dimensionless quan¬ 
tities. Treating all stocks Si as evolving under identical 
conditions is a strong assumption of the model but allows 
us to solely focus on the interplay between network and 
stock dynamics and its dependence on a few key param¬ 
eters. 

The effort is a time-dependent quantity assigned to 
each node Vi which defines its current behavioral pattern 
and changes through imitation of other nodes. On the 
one hand, nodes can adopt a high effort level > 1, 
causing each stock to converge to a stable fixed point 
5+ = 0 implying full depletion of the resource. Alter¬ 
natively, nodes can choose a low effort level E_ G (0,1) 
providing less harvest per unit time initially but avoid¬ 
ing depletion of the resource stocks since each individual 
stock Si then converges to a stable positive fixed point 
S- = 1 — E_ > 0. The two possible choices of effort level, 
E_ (low) and E_^ (high), are the same for all nodes and 
are parameterized by AE G (0,1) such that E_ = 1 — AE 
and = 1 + AE. At each time t there are N_ (t) nodes 
with Ei (t) = E_ and N_^ (t) = N — N_ (t) nodes with 
Ei(t) = E^. The effort then yields for each node Vi an 
individual harvest hi{t) = Si{t)Ei{t)^ which constitutes 
the second term in Eqn. Q. Erom now on we omit the 
explicit time dependence of the stocks 5^, efforts the 
adjacency matrix A and the number of low and high ef¬ 
fort nodes , in our notation. 

Initially, for each node Vi, an individual waiting time Ti 
is drawn at random from a Poissonian distribution with 
density 


p(T,)=T-iexp(-Ti/T), (2) 

which is a typical choice for modeling interaction rates in 
social systems [30] . T denotes the expected waiting time 
between two interactions initiated by the same node Vi. 
Starting from this: 

(i) The system as given in Eqn. 0 is integrated for¬ 
ward in time for the minimum of all current waiting 
times Ti. Then, for the corresponding node Vi (with 
the smallest Ti) a neighboring node Vj is drawn uni¬ 
formly at random. 

(ii) If the efforts Ei and Ej of Vi and Vj differ: 

(a) With a rewiring probability 0 < 0 < 1, breaks 
its link with Vj such that Aij = 1 becomes 
Aij = 0. Then, a new link between Vi and an¬ 
other randomly drawn node Vk with the same 
effort level {Ei = E^) is established such that 
Aik = 0 becomes Aik = 1- This network adap¬ 
tation process mimics generally observed tenden¬ 
cies to form clusters of individuals with similar 
behavior or social traits. Note that, in contrast 
to earlier work, rewiring only takes place if a ran¬ 
domly drawn neighbor Vj of Vi shows a different 
effort, e.g., behavioral pattern nni. 
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(b) If Vi does not adapt its neighborhood, imitation 
may happen instead (with probability 1—0). The 
difference in current harvest Ahij = hj — hi is 
computed and the node Vi imitates the current 
effort level of Vj with a probability given by a 
sigmoidal function p{Ei Ej) = p{Ahij) which 
generally is required to be monotonic and contin¬ 
uously differentiable. Additionally it must fulfill 
p{Ahij) 0 for Ahij —oc, p{Ahij) 1 
for Ahij CO and p(0) = 0.5. This repre¬ 
sents the increasing likelihood of imitation pro¬ 
cesses to take place with an increase in the ex¬ 
pected payoff difference m- For our model we 
set p{Ei Ej) = 0.5(tanh Ahij + 1) which obeys 
all of the above requirements. 

(iii) A new waiting time Ti is drawn at random for Vi 
according to Eq. and step (i) is repeated as long 
as the model has not reached a steady state. 

(iv) The model reaches (with probability one) a steady 
state at some time tf when the network divides into 
one or more components in each of which only one 
choice of effort level is left. 

Initially the two possible effort levels are distributed 
evenly among the nodes with ratios n_ (0) = N_{0)/N = 
n^(0) = N^{0)/N = 0.5. Initial stocks are set to 5^(0) = 
1 for all i = 1,..., A. In the following, we consider ini¬ 
tially Erdos-Renyi random networks with N = 400 nodes 
and a linking probability of p = k/{N — 1), where k = 20 
is the average degree of nodes in the network. 


III. STATIC NETWORK 

We first study the case of a static network structure 
with 0 = 0 (hence, modeling step (ii)(a) is not imple¬ 
mented at first) and simulate the model numerically for 
different combinations of T and AE. Erom this, we de¬ 
rive a macroscopic approximation of the model consti¬ 
tuted from a set of three coupled differential equations 
and show its good agreement with the numerical results. 


A. Numerical simulations 

Numerical simulations for different combinations of T 
and AE provide insights into this system’s dynamics. 
Eigure (A) shows the mean fraction f-{tf) of model 
runs that converge to a state where all nodes show a low 
effort Ei(tf) = E_ V i = 1,...,A (using an ensemble 
of n = 500 simulations). Eor small T (fast interactions) 
there is a high probability for the system to converge to 
a state where only nodes with a high effort level E^ are 
present. In this case all resource stocks converge to the 
stable fixed point 5+ = 0 and become fully depleted. 
With increasing T, the system’s expected equilibrium 
state undergoes a phase transition in /_(tj). Eor suf¬ 
ficiently large T (slow interactions), the system is likely 


to converge to a state where all nodes adopt the effort 
level E_ and all stocks converge to a stable fixed point 
5_ = 1 — £’>0. This indicates that the rate of interac¬ 
tions between nodes plays a crucial role in determining 
the system’s expected equilibrium state. 

The resulting dynamics can be qualitatively under¬ 
stood by considering the limiting cases of T ^ 0 and 
T ^ oc. In the first case, interactions between nodes 
are expected to happen very fast. Given that initially all 
stocks carry the same value 5^(0) = sq we expect that for 
t 1 the harvest h_ (h^) of nodes with low (high) effort 
follows h_{t ^ 1) oc E_so {h^{t ^ 1) oc E^sq). This im¬ 
plies that the difference in harvest between the two differ¬ 
ent types of nodes is expected as h_^ {t 1) — h_ {t 1) (x 
{E_^ — E_)so = 2AEso. If interactions happen very fast, 
the system likely converges into its equilibrium state at 
tf <^1. Since in this situation we expect h_^ > , nodes 

with low effort are more likely to imitate the high effort 
rather than the other way around and, hence, we expect 
^ 0 for T ^0 (as can be seen in Eig. 111(A)). 

In contrast, for T ^ oo we expect updates between 
nodes to happen preferably at times t ^ 1. In this case, 
the stocks of nodes with high (low) effort can be assumed 
to have already converged to a fixed point of 5+ = 0 
{s- = 1 — E_^ = AE) as interactions between nodes start 
to take place. Hence, the difference in harvest is expected 
as h_(t ^ 1) — h^(t ^ 1) = AE — A£^. Thus, for 
all AE G (0,1) the harvest of low-effort nodes exceeds 
that of nodes with high effort and the system is likely to 
converge into a state where all nodes show the low effort 
and, hence, f-{tf) 1 (red/dark area in Eig. (A) for 
high values of T). 

We note, that h_{t ^ 1) — hj^{t ^ 1) = AE — AE‘^ 
varies with AE. Specifically, in the limiting cases AE = 
0 and AE = 1 we find that the difference h_{t ^ 1) — 
h^(t ^ 1) = 0 vanishes and, hence, the system becomes 
equally likely to converge into either a state with only 
low-effort nodes or only high effort nodes present (see 
lower right corner and the shift of the transition point 
towards higher T with increasing AE in Eig. 0(A)). 


B. Macroscopic approximation 

Abstracting from pairwise microscopic interactions, we 
now look at the system from a macroscopic point of view. 
Assuming the network to be large and fully connected 
at first, the time evolution of the system’s state can be 
characterized by rate equations for three quantities: (1) 
the fraction of nodes n_ with effort level £_, (2) the mean 
resource stock p_ = {si\Ei = E_)i of nodes with effort 
level E_ and (3) the mean resource stock = {si\Ei = 
E^)i of nodes with effort level E^. The fraction of nodes 
with effort level E^ follows from = 1 — n_. 

The time evolution of n_ is governed by nodes that 
change from the low to the high effort level and vice versa. 
In particular, in the time interval (t, t -1- dt) an infinitesi¬ 
mal fraction of dn _ {dn_^^_) nodes change their effort 
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FIG. 1. (Color online) (A) The mean fraction /-(t/) of nu¬ 
merical simulations that converge to a state where all nodes 
show a low effort level Ei{tf) — E_ V z = 1,..., A" computed 
over n — 500 runs for different choices of T and AA for a 
static network with 0 = 0. (B) The value n_^ of the sta¬ 
ble fixed point for the fraction n_ of nodes with effort level 
E_ computed from Eqs. |T^-( 18). The dashed line indicates 
the critical expected waiting time Tc which separates the two 
regimes (predominance of nodes using F0 (yellow/light) and 
E_ (red/dark)). 


from E_ {E^) to E^ {E_) which decreases (increases) the 
fraction of nodes with low effort n _, 

dn_ = dn_^^_ — dn _(3) 

The interactions between nodes that govern the rates of 
changes in effort are driven by the following quantities: 

1. The expected waiting time T for a node Vi to inter¬ 
act with a randomly drawn neighboring node Vj. 
Correspondingly, the rate of node interactions is 
taken to be r = 1/T. 


Ej = E^ {Ej = E_). This quantity is then depen¬ 
dent on the expected stocks at nodes with low and 
high effort, which is derived below in detail. 

This yields dn _and dn^^_ as the product of all three 

factors introduced above, 

dn = n_rn^p ^^dt (4) 

dn^^_ = n_^rn_p_^^_dt ( 5 ) 

dn_ 

^ ~di -P-^+)' W 

The two quantities still remaining to be evaluated, are 

the expected probabilities p_^^_ {p _,+ ) for nodes with a 

high (low) effort level to change to the opposite level. It 
is obtained as the expected probability for nodes in the 
network to take up its neighbor’s effort, 

P+^- = {P{^j ^k) I Ej = E^, Ek = E_)j^k 

= 0.5(tanh(A/zj/c | Ej = E^ = E_))j^k + 0-5 
= 0.5(A/z^-/c I Ej = E^, Ek = E_)j^k +0.5 
= 0.5(-E_ [sk I Ek = E_)k — E_^{sj I Ej = E^)j) 

+ 0.5 

= 0.5(£L/i_ -A^/i^) + 0.5 (7) 

p_^^ =0.b{E^p_^ -E_p_) + 0.5. (8) 

Here, we performed a linear expansion of the hyperbolic 
tangent, tanhx = x + O(x^), assuming that differences 
in harvest remain small. 

The time evolution of either of the two average stocks 
p_ and is governed by two terms. First, each individ¬ 
ual stock Si follows the logistic growth model and so do 
the average quantities. Second, the value of each of the 
two average stocks changes according to the fact that the 
nodes modify their effort from E_ to E^ and vice versa 
during the time interval (t, t -\- dt). This yields 


2. If a node Vi interacts with its neighboring node Vj, 
an imitation of effort only takes place if Ei ^ Ej. 
Hence, for a node vi with Ei = E_ {Ei = E^) there 
is to define a probability {P^ ) that a randomly 
drawn neighboring node Vj has Ej = E^ {Ej = E_). 
Since a large fully connected network is assumed, 
this probability is given exactly by the current frac¬ 
tion (n_) of nodes with high (low) effort E^ {E _) 

and, hence, P^ = n_^ {P^ = ti_). 

3. If a node Vi with Ei = E_ {Ei = E_^) interacts with 
a neighboring node Vj with Ej = E^ {Ej = A ), 

there is a probability p _that Vi takes up 

the effort level Ej of Vj. This probability is gov¬ 
erned by the difference in harvest Ahij between 
Vj and Vi. For the macroscopic description, the 
individual pairwise interactions are replaced by ag¬ 
gregated quantities. Therefore p_(p^^_) is com¬ 

puted as the expected probability for a node Vi with 
low (high) effort to adopt the high (low) effort given 
that it interacts with a node Vj that currently has 


dp_ — d{sk I Ek — E_'j k 

— {dsk I Ek — E_^k 

— dt {sk{^ Sk^ EkSk I Ek — E_^k + 

= dtp_ — dt{s\ \ Ek = E_)k — dtE_p_ + S_ 

= dt{p_ {1 — p_ — E_) — + 5_ (9) 

- EV - +^) + <5+. (10) 

Here p^^ and p^^ denote the variances in the two types 
of stocks. 5_ ((5^) indicate the net change in the average 
stock as nodes with high (low) effort change their effort 
to the opposite choice during (t,t + dt). The fraction of 

nodes dn_^^_ {dn _^^) that change their effort from E_^ 

to E_ {E_ to E^) during (t, t + dt) is assumed to be small 
compared to the fraction of nodes which already hold the 

low (high) effort, dn^^_ n_ {dn _<C n^). Hence, 

the respective contribution to the dynamics of p_ {p^) 
as nodes change their effort is also assumed to be small, 

dn_^^_p_^ n_p_ {dn _ ^^p_ This allows for a 

first-order expansion of the stock’s time evolution, such 
that 
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S_ — 


(n_ — dn _ 

n_ — dn _+ dn^^_ 

{n_ — dn _+ dn_^^_jii_^ 


+ 


n_ dn_^^ + dn^^_ (dn_^^ ,dn^^_ )=(o,o) 

— {n_ — dn _+ dn_^^_) + ((n_ — dn _+ dn_^^_fi_^) 


{n_ — dn _+ dn^^_)‘^ 

/i^(n_ — dn _+ dn_^^_) — {{n_ — dn _,+ )//_ + dn_^^_jii_^) 


dn 


(dn ^ , ,dn. ) = (0,0) 


(n_ — dn _+ dn_^^_)‘^ 


dn^ 


(dn ^ , ,dn, ) = (0,0) 


n_ii_ —ld_n_ -\- n_jd_ ^ —Ti_ii_ M+ ~ /^_ 


+ ■ 

n_ n^ 

= (aV — ld_)n_^Tp_^^_dt 
= (/^_ — ld_^)'n_rp _ ^_^dt. 


-dn _+ 


n^ 


-dn _ = /i_ + 


dn^ 


n 


(11) 

( 12 ) 

( 13 ) 


Putting this back into @ and ( P!q| ) yields 

dp_ = dt{p_ (1 — /i_ — E_) — 

+ dt{p^ - p_ )n^rp^^_ ( 14 ) 

dfi^ = dt{p^{l - E^) - 

Edt{p_ - p^)n_Tp_^^. ( 15 ) 

In the scope of this work, in to order to close the set 
of equations that describe the systems dynamics, we as¬ 
sume the respective variances and to vanish. 
Taking into account higher moments in the dynamics of 
the stocks and investigate its influence on the resulting 
fixed points remains as a task for future research. 

In summary, we find a set of three coupled ordinary 


differential equations that define the time evolution of 
the static network model: 


— = rn+n_ (16) 

dP- / X / X / X 

(1 - M- - -B-) + 'r(M+ - )n+P+^- (17) 

^ = li+(l - M+ - ^1) + t { ii _ - Ai+)n_p_^+. (18) 


C. Fixed points and stability 


We obtain all fixed points Pi = {n_^p_ p ^of the 
dynamical system given in Eqs. (16)-(18) as: 


Pi = 
P2 = 

P3 = 

Pa = 

P5 = 


^ 1-E_- 0.5r 

4-0 = °’ ^-0= l + 0.5rE ’ ^-0 = 0 

^ l-E.- 0 . 5 t 

n_„ = l, ^_„ = 0, 1 ^. 0 = i + o.5rK 


2{E_ 


1-0.5r 


E^-l) 


E_+E ' ■^+ 1 — 0 . 5 r ^ 1 — 0 . 5 t 

0 “ JeI t; ’ l^-o~ E + E ’ ^+o~^-' 


^(fr-1) 


E_+E^ 


^-o-r ii-o-i--Ei> i^+o-^ + v(^) +- 


2a 


n_o = 1, A^-o = 1 - £1, 1^+0 = ^ - \/ ( ^ ) + - 


2a 


a = 0.5{-2 - E^t) 

b = l-E^+ 0 . 5 t ((1 - E_)E^ +E_- E^ - 1 ) 


c = 0.5 t(1 - E_ ){E_ - E^ - 1). 


( 19 ) 

( 20 ) 

( 21 ) 

( 22 ) 

( 23 ) 
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In addition, there exists a manifold which also satisfies 
^ ^ ^ = 0 and is defined by 

Pa = (n_o = a, I 1 _^ = 0 , = 0 ), a e [ 0 , 1 ] (24) 

For all fixed points given above we compute the largest 
eigenvalue Ai of the corresponding Jacobian matrix eval¬ 
uated at the respective point. Only the two fixed points 


P 3 and P 4 have a negative largest eigenvalue Ai < 0 
and, hence, are stable for choices of parameters AP and 
T > 0.5 (note that again: P_ = 1 — AP, P^ = 1 + AP 
and r = 1/T) (Fig. [^. 

To investigate the system’s dynamics in the regime T < 
0.5, the stability on the 1-dimensional manifold defined 
by all points that fulfill Eq. (24) is assessed. Analytically 
computing the three eigenvalues of the Jacobian matrix 
on the manifold as a function of the parameter a yields 


Ao = 0 


(25) 


\±{a) = 1 - ^ ± \\j2aE^T - 2aE_T ■ 


■ E^ - 2E^E_ - E^t - 


-E^ 


■ E T + —. 
4 


(26) 


A first observation is that A+((a) > A_((a) holds. Since 
Aq = 0 , it is obvious that not all eigenvalues can be neg¬ 
ative. However, if Aq = 0 is the largest eigenvalue of the 
system, all choices of a for which A+((a) < Aq define a 
center manifold. 


Hence, 


A+(<a) < 0 if (a < - — TAP 


iy(a) = (nso = a, fiso = 0, fino = 0) 


OL G 


0, - - TAE 


(27) 

(28) 


defines a center manifold where the system’s stability 
cannot be assessed by linear stability analysis. A de¬ 
tailed study of the system’s stability in this regime is be¬ 
yond the scope of this work and not necessarily needed to 
understand the general dynamics of the macroscopic de¬ 
scription proposed here. Numerically integrating the sys¬ 
tem for choices of parameters taken from the center man¬ 
ifold, however, reveals good agreement between the mi¬ 
croscopic and macroscopic model representation (Fig.[^. 
An investigation by means of a higher-order stability 
analysis might, however, yield further insights into the 
processes that cause both resource stocks /i_Q = /i+Q = 0 
to be fully depleted in the regime of the center manifold. 

In conclusion, we note that for each choice of param¬ 
eters only one of the fixed points P 3 and P 4 can be the 
unique stable fixed point of the system (Fig. §. Fig- 
ure (B) displays the value of the stable fixed point’s 
n_ Q-component as a function of T and AP. The results 
are in good agreement with the numerical findings (Fig.[^ 
(A)). Due to the first-order approximation, the transition 
from a predominance of nodes with P^ to nodes with P_ 
with increasing T is not as sharp as for the numerical 
simulations. However, a good estimate for the critical 
value Tc of T at which the transition takes place can be 
found by setting n_^{Tc) = 0.5 in Eq. (21) which yields 

Tc{AE) = 2 -^e‘^ (dashed line in Eig. [l). 



FIG. 2. (Color online) The largest eigenvalue Ai for the two 
fixed points P 3 (A) and P 4 (C) (see also Eqs. ( [^ and (22)) 
depending on AP and T. The black area in (B) indicates 
the domain in parameter space for which Ai computed for 
P 3 is negative and, hence, P 3 is stable. (D) shows the same 
properties for P 4 . The regimes for which either of the two 
fixed points is stable are complementary. Further it should 
be noted that for T <0.5 neither of the two fixed points is 
stable, but center manifold as given in Eqn. (28) exists in this 
regime. 


IV. ADAPTIVE NETWORK 


In the following, we consider additionally network 
adaptation processes with 0 > 0 (hence, modeling steps 
(ii)(a) and (b) both take place with a relative frequency 
depending on the rewiring probability 0). For all results 
presented from here on, the two available choices of effort 
levels are fixed by setting AP = 0.5. 
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A. Numerical simulations 

Numerical simulations with the same initial conditions 
as in the static case for different combinations of (j) and 
T reveal a division of the parameter space into regimes 
of different expected outcomes as the model reaches its 
steady state (Fig. [^(A)). As for 0 = 0, fast interactions 
(i.e., low values of T) lead to a large fraction of nodes 
carrying . The transition between the two behavioral 
patterns with increasing T remains sharp. However, de¬ 
pending on the choice of 0, the value of the critical wait¬ 
ing time Tc, at which the system transfers from a state 
with a predominance of nodes with low effort to a state 
with a predominance of nodes with high effort, decreases 
with increasing (j). Conversely, this implies that for all 
T > 0.3 there is an appropriate choice of 0 G [0ci:0c2] 
so that all nodes are likely to adopt the effort level E _. 
In the limiting case of 0 = 1 the expected fraction of 
nodes with E_ equals the initial fraction n_ (0) = 0.5 for 
all choices of T due to the network’s fragmentation into 
components of nodes sharing the same effort. 


B. Macroscopic approximation 


The macroscopic approximations (16)-(18) can be ex¬ 


tended to also include the effects of network rewiring. For 
this, we introduce two additional variables describing the 
macroscopic state of the network. The time evolution of 
the fraction of nodes n_ with low effort is recalled (anal¬ 
ogously to Eq. ®) as 

^ = T{n^P^p+^_ - n_PtP-^+)- (29) 


Given that a node Vi initializes an interaction and the 
randomly drawn neighboring node Vj employs a different 
effort, Ei 7 ^ Ej, there exists the adaptive rewiring prob¬ 
ability (j) G [0,1] for Vi to break its connection with Vj 
and establish a link with another randomly drawn node 
Vk in the network that is employing the same effort as 
Vi {Ek = Ei) and is not yet connected to node Vi. With 
probability 1 — 0, imitation of efforts takes place which 
has already been implemented in the macroscopic de¬ 
scription of the static network. To account for the adap¬ 
tive rewiring process, the interaction rate r needs to be 
refined such that it no longer represents the rate of node 
interactions alone, but the rate of interactions which lead 
to imitation. 



FIG. 3. (Color online) (A) Mean fraction of nodes f-{tf) 
with effort level = 1 — AE = 0.5 for different choices of 
T and 0 obtained from an ensemble of n = 500 numerical 
simulations as the system reaches its steady state. (B) Value 
of the stable fixed point for the fraction of nodes with effort 
level E_ computed from the set of differential equations (j^- 


For adaptive rewiring to take place, the network cannot 
be fully connected. Therefore, the previous definitions 
of = n_ and P_ = n_^ for two nodes of different ef¬ 
fort to interact no longer hold for the derivations to be 
performed here. 

The total number of M links in the network splits into 
M_ {M_^) links connecting two nodes with low (high) ef¬ 
fort and M_^_ links connecting two nodes of different ef¬ 
forts, such that 


M = ^ = M_ 

dM dM_ dM^ dM^_ 

dt dt dt dt 


(32) 

(33) 


Additionally let 


K 


2M_ 

~jr 


(34) 


denote for nodes with low effort the average number of 
neighbors with the same effort. Likewise, 


K 



(35) 


r = 


T 


(30) 


Likewise the ratio p of all node interactions that lead to 
adaptive rewiring needs to be defined. Since each node 
is expected to interact at a rate 1/T it follows that 


represents for nodes with low effort the average number 
of neighbors with high effort. These two quantities con¬ 
stitute the average degree of nodes with low effort as 


K_=K_+ iV 


+ 2M_ 

V 


(36) 


P = 


T' 


(31) 


Likewise the average degree of nodes with high effort 
























is obtained from 


+ 

2M^ 


N 


(37) 


Af 


K~ 

+ 

H— 

(38) 


M^ + 2M^ 



K ■ 

(39) 


For a node Vi currently having a low effort Ei = E_ the 

probability (vi) to draw a neighbor vj with different 
effort at random is given as 




k^(vi) 

k(vi) ' 


(40) 


(vi) is the number of neighbors of node Vi that employ 
the high effort. k(vi) denotes the degree of node Vi. Since 
for the macroscopic description the pairwise microscopic 
interactions between nodes are approximated by the av¬ 
erage dynamics, we compute the average probability P 
for a node Vi with low effort to interact with a node em¬ 
ploying the high effort. Since the network is initialized as 
an Erdos-Renyi random network and it is further equally 
likely for all nodes with the same effort to connect to 
or disconnect from other nodes by random rewiring, we 
perform a heterogeneous mean-field approximation and 
assume the average degree k{vi) to be the same for all 
nodes with low effort, k{vi) =K V i G {1 ,..., | = 

£L} [3ll[32]. Thus 



FIG. 4. (Color online) Illustration of the influence of the 
imitation of effort on the different numbers of link types in 
the network. A node Vi with the high effort Ei = E_^ (in¬ 
dicated in orange) interacts with a node Vj with low effort 
Ej = E_ (red/dark). Node Vi may then imitate the effort of 
node Vj, Ei ^ E_ . The number of links between nodes with 
low (high) effort M_ (M^) then increases (decreases) by the 

number k_^{vi) (/q(r’i)) of neighbors of Vi that show the low 
(high) effort. 


and is fully determined by the per node densities of links 
and m _. 

Generally, the time evolution of the total number of 
links between nodes of low effort is governed by imita¬ 
tion and adaptation. First, we focus on the process of 
adaptation. Since links between nodes of the same effort 
can only be established but not removed via the process 
of adaptation, the contribution of this process to the total 
number of links between low-effort nodes M_ only causes 
it to increase. This positive contribution is 


dM_ 

dt 


pN_P'_ 


(45) 


{P_ {vi) \Ei=E_)i = 


' k_{vi) 

K_ 

M^_ 


Ei=E 


2M_ + M.. 


k^{vi) 

k{vi) 



E 

HI 


(41) 


Instead of the actual number of M links in the network 
we define the corresponding per node link density 


M 

^=N = 


N 


M_ iir 


2 


■ m 


■ m. 


(42) 


which is independent of the number of nodes N in the 
network. The probability for a node with low (high) ef¬ 
fort to interact with a node of high (low) effort is then 
given by 




2m_ + 
m^_ 
2m^ + 


(43) 

(44) 


and is explained as follows: In each time interval (t, tpdt) 
there is a total number of N_ nodes, which with prob¬ 
ability p initiate an interaction that leads to adaptive 
rewiring. Adaptive rewiring then takes place if a ran¬ 
domly drawn neighbor Vj of the considered node Vi em¬ 
ploys the high effort. As defined above, this happens with 
probability P . 

The second contribution to the time evolution of M_ 
is given by imitation, which takes place at rate r. Gen¬ 
erally, there is one term causing an increase in links be¬ 
tween nodes with low effort and one term causing its 
decrease. First, assume a node vj with Ej = E_^ to imi¬ 
tate the low effort E_ from one of its neighboring nodes 
Vi with Ei = E_. The number of links between nodes 
of low effort then increases by the number k^ (vj) of all 
neighbors of node Vj that employ the low effort (Fig. [^. 
Again, by performing a heterogeneous mean-field approx¬ 
imation and assuming the number of neighbors for indi¬ 
vidual nodes to be represented by the respective average 
number of neighbors, we set 

K{vj)=K = ^- ( 46 ) 

Now, it holds that each of the N_^ nodes with high effort 
interacts with a node of low effort with probability P^ 
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at rate r. Then, with probability p_^^_ a node with high 
effort takes up the low effort. This causes the number of 
links between pairs of nodes with low effort to increase by 
the number of neighbors with low effort of the formerly 
high-effort node, 


(47) 

A third term that governs the time evolution of M_ 
is given by its decrease caused by nodes with low effort 
that imitate the high effort. If a node Vi with the low 
effort Ei = E_ interacts with a node Vj having the high 
effort Ej = E_^ and Vi then imitates the effort of Vj^ the 
total number of links connecting two nodes with low ef¬ 
fort decreases by the number of Vi^s neighbors Vk that 
are showing the low effort Ej^ = E_ as well. Following 
from an analogous argument as given above, this num¬ 
ber is given by k_ (vi). Again we assume the number of 
neighbors Vk with E^ = E_ of a node Vi with Ei = E_ to 
be approximated by its average. 


k_{vj) = K_ 


2M 


(48) 


With rate r each of the N_ nodes with low effort interacts 
with a node showing the high effort E_^ with probability 

P\ With probability p _a node with low effort imi¬ 

tates the high effort which causes a decrease in M_ by the 

average number of low-effort neighbors K of the node 
that is imitating the high effort. 


dM_ 

dt 


-tN_P^P_ 


K 


(49) 


Putting together Eqs. (45), m) and (49) gives the time 


evolution of the number of lin 
effort as 


^s between nodes of low 


dM_ 

dt 


r(7V+P>+^_< -iV 

+ pN_p\ (50) 


Plugging the definitions of K (Eq. (HI) and 
(Eq. ([ 3 ^) into Eqn. (50) and normalizing with the to¬ 
tal number of nodes N yields the time evolution of the 
per node density of links between nodes of low effort 


dm_ 

dt 


T{P^p^^_m_^_ 


2P^p _ ^^m_) pn_P^ (51) 


which is again independent of N. Due to the symmetry of 
the system, the time evolution of the per node density m_^ 
of links between nodes with high effort then immediately 
follows as 


dm^ 

dt 


- 2P^ P+^_m+) + pn^P^ . (52) 


Eor the time evolution of the average stock of nodes 
with low and high effort p_ and we already found in 
Eqs. and ( [1Q| ) that 


dp_ = dt{p_ (1 — /i_ — E_) — + 6_ (53) 

- E_^) - pP ) + < 5 ^. ( 54 ) 

The general forms of S_ and are (see Eq. ( p!^ and ( p!^ ) 


= 




p^ - p 


dn, 


-dn 


(55) 

(56) 


Eor the case of an adaptive network, dn^^_ {dn _^^) is 

given by the first (second) term in Eq. (29): 


<5_ = ^~ Tn^Pp+^. (57) 

6, = — ^Tn_P^p , (58) 

n_^ 


with the probabilities P^ and P^ (Eqs. ( [4^ and (44)) as 
defined above and p_^^_ and P_^ being the same as for 
the static model (Eqs. 0 and ^). 

To summarize, the set of five coupled differential equa¬ 
tions that represent the adaptive network model’s macro¬ 
scopic dynamics is given as 


= T(n^ P^ p^^_ - n_ P* p_^^) (59) 

—^=r{P^p_^^_m_^_—2P_p _ P pn_P_ (60) 

dm, , + - , - , , 

— = t{P_ p_^^m^_ - 2P^ P+^_TO+) + pn^P^ (61) 

^ - 1^- - E_) + - p_)P^ p^^_ (62) 

^ = Ai+(1 - M+ - -B+) +'r^(/i_ - (63) 


It is important to note that in most previous works on 
adaptive networks a closed set of macroscopic equations 
is obtained by assuming that links in the network are 
drawn at random and interactions take place between 
nodes that are connected by them [Hiisa]- In this work 
nodes, not links, are randomly drawn and initiate an in¬ 
teraction with its neighboring nodes. This subtle dif¬ 
ference changes the effective time scale of the system. 
Specifically, in our model only a maximum of N out of 
all M links are affected by interactions between nodes 
during the same time as all M links would be considered 
if interactions take place by randomly drawing links in 
the network. In other words, in our model it takes M/N 
times longer to achieve the same number of updates, as 
one would obtain by considering per-link interactions. 

Eor this system, the stable fixed point’s n_ ^-component 
can be obtained numerically for different combinations 




















10 


of <p and T (Fig. (B)). The results are again in good 
agreement with the numerical simulations and imply that 
for every choice of T > 0 there actually exists an ap¬ 
propriate choice of 0 G [ 0 ci 5 0 c 2 ] so that all nodes are 
likely to adopt the effort level E_. The lower bound 
of the optimal rewiring probability can be obtained 
by utilizing Eq. (21) and the linear relationship between 
0ci and T for a fixed rate of social updates r that 
lead to imitation as given in Eqn. (30). We thus find 

2-2AE‘^ 


2AE‘^ 


(l^eATAE) = V 0 < T < ^ 

and (^d {T, AE) = 0 otherwise. The upper bound (/)c 2 at 
which the network fragments is obtained from a numer¬ 
ical bifurcation analysis as 0C2 ~ 0.89. The result is in 
good agreement with previous findings on the fragmen¬ 
tation threshold in adaptive networks for similar average 
degree k [34l[35]. We find, however, that the computed 
fragmentation threshold 0 C 2 is larger than what is ex¬ 
pected from the numerical simulations (Eig.[^(A)). This 
can either be due to the fact that moment closure as 
well as mean-field approximations are known to provide 
only rough estimates of the fragmentation threshold [33] 
or because finite size effects in the numerical simulations 
cause the system to fragment for smaller values of 0 than 
it would be expected for the limiting case N ^ oo that is 
considered in the macroscopic approximations. A more 
detailed study of the network fragmentation and the cor¬ 
responding threshold 0 C 2 is a subject of future research. 


C. Consistency between approximations 


To illustrate the consistency of the set of differential 
equations describing the static setting ([l^-( 18) and the 


adaptive case (59)-(63), we set 0 = 0 in the latter, com¬ 


pute its fixed points numerically and compare them with 
the static setting’s fixed points (21) and (22). Eigurej^ 
(A-C) shows the different components of the stable fixed 
points as a function of the control parameter T for a 
fixed AE = 0.5. The components n_Q, and /i_Q align 
perfectly well for the static and the adaptive case. The 
gray shaded area in Eig. (A) indicates a center mani¬ 
fold for which the system’s stability cannot be assessed 
by standard linear stability analysis. However, numer¬ 
ically integrating the set of differential equations yields 
the expected behavior of n_ (0) ^ 0 as T ^0. Eig- 
ure[^(D) displays again the n_ q- component of the adap¬ 
tive model’s stable fixed point for 0 = 0 and different 
combinations of T and AE. The results match those of 
Eig.[^(B). Hence, the system of dynamic equations (59)- 
( [^ can be interpreted as a consistent generalization of 
Eqs. ([^-(@. 


V. CONCLUSIONS 

We have introduced a model to describe emerging 
structure formation from the interplay of dynamics of 




,1.00 


0.75 


0.50 


0.25 


0.00 


FIG. 5. (Color online) (A)-(C) The dependence of the adap¬ 
tive (solid lines) and static model’s (transparent scatter) sta¬ 
ble fixed point on the expected waiting time T for fixed pa¬ 
rameters 0 = 0 and AE — 0.5. (D) The adaptive model’s 
stable fixed point’s n_ ^-component indicating the fraction of 
nodes with effort E_ in the consensus state as a function of 
the two parameters T and AE again for 0 = 0. The dashed 
line indicates the value of the critical wait in g tim e E obtained 
from the set of differential equations ( 16)-([T^. 


and on networks manifested by the co-evolution of so¬ 
cial dynamics on the one hand and resource dynamics 
on the other hand. An adaptive voter model has been 
coupled to a set of logistic growth models, such that the 
state of the dynamic variables influences the imitation 
(i.e. social trait adoption) processes in the underlying so¬ 
cial network which take place according to differences in 
harvest or payoff. We have derived rate equations for the 
system’s macroscopic variables and demonstrated that 
the resulting system of differential equations yields sta¬ 
ble fixed points which are in good agreement with the 
results from numerical simulations. 


Our paradigmatic example illustrates that the inter¬ 
play between both types of network dynamics gives rise 
to a variety of new phenomena, which have not been 
observed so far when only studying either of the two as¬ 
pects. We have mainly found that the rate of interactions 
in the network determines the expected linear stability of 
the growth model’s fixed points. However, for each choice 
of interaction rate there exists an appropriate range of the 
adaptive rewiring frequency so that the expected fraction 
of, e.g., nodes with effort E_ can be maximized. Notably, 
the subset of differential equations (59)-(61) provides a 
general description of imitation and adaptation dynamics 
on a social network with binary states of nodes and sym¬ 
metric imitation rules. Hence, it is applicable to study 
many other problems as long as the imitation probabil¬ 
ities p _and , which do not have to be constant 

for all times, are chosen appropriately. 

The proposed model also raises questions that need 
to be addressed in future research. In the course of 
the macroscopic approximation we have assumed all mo¬ 
ments of higher order in stocks and network structure to 
vanish such that the set of differential equations could 
be closed. The results have been shown to be in good 
agreement with numerical simulations. However, a more 
in-depth analysis of whether the inclusion of higher order 
moments would enable us to reproduce the steep transi- 
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tion between the two regimes of predominance of low- or 
high-effort nodes remains a relevant research questions. 
We also aim to estimate more thoroughly the critical 
waiting time Tc at which the observed phase transition 
takes place and therefore investigate the expected time 
at which the low effort provides more harvest than the 
high effort given that no interaction between the nodes 
took place so far. Finally, we aim to obtain data from 
agricultural studies on, e.g., water usage or harvest ex¬ 
ploitation of resources, to test the findings and insights 
that we have obtained from our co-evolutionary model 
with respect to real-world phenomena. 
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